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ABSTRACT 

A  domain  decomposition  algorithm  suitable  for  the  efficient  and  accurate  solu¬ 
tion  of  a  parabolic  reaction  convection  diffusion  equation  with  small  parameter  on 
the  diffusion  term  is  presented.  Convergence  is  established  via  maximum  principle 
arguments.  The  equation  arises  in  the  modeling  of  laminar  transonic  flow.  Decom¬ 
position  into  subdmomains  is  accomplished  via  singular  perturbation  analysis  which 
dictates  regions  where  certain  reduced  equations  may  be  solved  in  place  of  the  full 
equation,  effectively  preconditioning  the  problem.  This  paper  concentrates  on  the  the 
theoretical  basis  of  the  method,  establishing  local  and  global  a  priori  error  bounds. 
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1.  Introduction.  In  this  paper  a  domain  decomposition  algorithm  for  the  solu¬ 
tion  of 

(1)  P[tt]  :=  U(  -I-  ttu,  -  ctt*,  —  ru  =  0, 

where  6  is  a  small  positive  parameter  is  presented.  This  equation  contains  many  of 
the  properties  that  make  the  gasdynamic  equations  difficult  to  solve;  namely,  it  is 
capable  of  modeling  rapid  variations  such  as  shocks  and  boundary  layers.  A  priori 
error  bounds  are  obtained  using  asymptotic  analysis,  and  are  verified  via  maximum 
principle  arguments.  The  analysis  also  identifies  parallelism  intrinsic  in  the  physics 
of  the  problem.  This  parallelism  may  be  exploited  by  the  particular  numerical  meth¬ 
ods,  allowing  efficient  use  of  parallel  arcUtectures.  This  paper  concentrates  on  the 
theoretical  basis  for  the  method,  discussions  of  the  numerical  details  may  be  found  in 
[10]  and  [15]. 

The  method  presented  here  is  appropriate  for  certain  problems  arising  when  mod¬ 
eling  laminar  transonic  fiow,  such  as  through  a  duct  of  variable  width.  When  modeling 
transonic  flow,  except  in  regions  of  rapid  variation  such  as  in  shocks  and  boundary 
layers,  convection  and/or  reaction  terms  dominate  over  diffusion.  The  reaction  term 
may,  for  example,  arise  from  the  effects  of  a  variable  cross  sectional  area  in  a  duct, 
thus  this  not  a  reacting  flow.  Asymptotic  analysis  identifies  the  regions  where  the 
solution  behaves  different,  subdividing  the  domain  into  the  following  two  types  of  re¬ 
gions:  regions  where  the  solution  is  smooth,  where  a  reduced  equation  may  be  solved; 
and  regions  of  rapid  variations,  such  as  in  a  neighborhood  of  a  shock,  where  the  full 
equation  must  be  solved.  The  domain  decomposition  is  independent  of  the  choice  of 
numerical  schemes  for  the  subdomains,  hence  the  numerics  will  be  discussed  in  this 
paper  only  briefly.  In  addition  to  dictating  the  domain  decomposition,  asymptotics 
also  provides  a  means  of  approximating  solutions  to  the  problems  in  the  subdomains. 
In  this  way,  a  set  of  simplified  problems  is  obtained  that  is  better  conditioned  for  nu¬ 
merical  computations.  The  domain  decomposition  and  preconditionings  are  reflected 
in  the  theorems  presented  herein. 

The  asymptotic  analysis  involves  the  derivation  of  analytic  upper  and  lower 
bounds  on  the  solution,  and  is  performed  in  the  style  of  Howes  [4,5,6].  The  method  is 
capable  of  obtaining  solutions  to  (1)  when  the  shock  is  not  stationary,  thus  extending 
Howes’  studies  [7,8]  into  the  time-dependent  repme. 

The  method  is  an  iterative  technique.  In  Section  3  the  domain  decomposition 
and  some  preconditionings  are  presented.  This  includes  an  error  analysis  of  the  pre¬ 
conditionings  and  the  theoretical  basis  of  the  domain  decomposition.  In  Section  5, 
the  method  is  summarized  by  outlining  the  algorithm. 

2.  The  Quasilinear  Problem.  Consider  the  behavior  of  the  solution  of  the 
quasilinear  parabolic  equation  (1)  on  the  domain 

(2) 

subject  to 

(3)  tt(x,0)  =  Tf(x),  0<x<6; 
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(4) 


tt(0, *)  =  a(t),  0<t<T\  and 


(5)  u(&,0  =  /9(*),  0<t<T. 

The  portion  of  the  boundary  along  which  the  data  is  specified  is  denoted  by 

n  :=  {(*,0|o  <  *  <  6,  t  =  o}f|{(*.0|o  <  t  <  r,  X  =  0,6>. 

For  the  sake  of  simplicity,  it  is  assumed  that  all  boundaries  are  inflow  boundaries, 
that  is,  o(t)  >  tto  >  0  and  0{t)  <  0o  <0. 

The  reaction  term  may,  for  example,  arise  from  the  effects  of  a  variable  cross 
sectional  area  in  a  duct.  Howes  [8]  discusses  the  case  when  r(z)  =  — a'(z)/a(z),  where 
a(z)  is  the  width  of  the  duct  (see  Figure  1).  The  coefficient  r  is  assumed  to  be  bounded 


A 


Fig.  1.  Variable  width  daet. 


with  bounded  derivatives. 

It  is  assumed  that  the  boundary  data  are  sufficiently  smooth  so  that  the  solution 
to  (1)  is  uniquely  defined  (for  example,  see  [2]).  We  are  interested  in  the  formation 
of  shocks,  thus  the  data  is  assumed  to  be  continuous.  For  example,  the  compatibility 
conditions 

(6)  a(0)  =  7(0),  and  7(6)  =  ^(0), 

must  be  satisfied.  In  addition,  it  is  assumed  that  the  first  derivatives  of  the  solution 
to  the  reduced  (e  =  0)  equation 

(7)  Po(t^l  :=  Ut  +  UU,-rU^0 

are  continuous  except  along  the  shock.  This  requires,  for  example, 

^  -  n  =  0,  for  (i,  t)  =  (0, 0); 


(8) 
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(9) 


^  +  *7^  -  fTf  =  0,  for  (x,  0  =  (6, 0). 


The  problem  is  essumed  to  be  nondimensionalixed  such  that  the  diffusion  co¬ 
efficient  c  is  inversely  proportional  to  the  Reynolds  number  (see  (13]).  Based  on 
free-stream  conditions  in  transonic  flow,  the  Reynolds  niunber  for  this  problem  is 
large.  Therefore,  it  is  appropriate  to  exploit  the  smallness  of  the  positive  parameter 
e  in  the  analysis. 

8.  Analysis  of  the  Precondltlonings.  The  method  exploits  preconditionings 
to  obtain  computational  efficiency  and  accuracy.  The  meaning  of  preconditioning  is 
broader  than  the  ustial  meaning  applied  in  the  linear  algebra  setting.  In  this  setting, 
a  problem  is  preconditioned  if  it  is  more  tractable  with  respect  to  numerical  compu¬ 
tations.  For  example,  in  this  section  a  preconditioner  based  on  a  physically  motivated 
domain  decomposition  is  discussed.  Regions  where  the  solution  behaves  differently 
are  identified.  The  problem  is  better  conditioned  because  each  numerical  method 
used  now  needs  only  capture  one  type  of  behavior  in  the  solution.  In  addition  to 
the  domain  decomposition,  a  preconditioner  involving  a  transformation  of  the  spatial 
coordinate  and  a  preconditioner  involving  a  modification  of  the  governing  equation 
will  be  presented. 

The  domain  decomposition  and  the  use  of  the  reduced  equation  are  closely  related. 
Asymptotic  analysis  identifies  two  types  of  regions.  In  the  outer  regions,  the  solution 
is  slowly  varying  and  the  eu„  term  is  small.  Thus,  in  the  outer  region  subdomains, 
the  governing  equation  is  modified  by  dropping  this  term  with  minor  effects  on  the 
error.  The  solution  to  the  reduced  equation  will  be  described  next. 

Let  17  be  a  weak  solution  of  (7)  with  boundary  data  (3-5),  which  is  a  solution 
to  (1)  in  the  limit  as  €  i  0.  For  the  analysis  here,  assume  that  U  has  a  single  shock. 
Let  the  path  of  the  shock  be  ipven  by  the  curve  (z,t)  =  (r(t),().  The  initial  and 
boimdary  data  are  assumed  to  be  smooth;  thus,  the  shock  does  not  exist  at  t  =  0. 
Rather,  F  is  assumed  to  be  imdefined  for  t  <  t**,  where  t  =  is  the  time  U  becomes 
discontinuous.  It  is  natural  to  describe  U  in  terms  of  the  following  functions: 

Uq{x,  t)  for  0  <  I  < 

I7(x,()  =  •  Ui{x,t)  for  X  <  T~^(i)  and  t>t^ 

UriXft)  for  X  >  r~^(f)  and  t  > 

For  analytic  methods  to  choose  F,  see  Whitham  [19]  or  Kevorkian  and  Cole  [9]. 

The  shocks  in  the  system  are  assumed  to  be  physical;  thus,  the  solution  will  satisfy 
the  entropy  condition 

(10)  t7,(F(«),0>*>ff,(F(f).*). 

where  the  speed  s  of  the  shock  is  given  by  the  Rankine-Hugoniot  jump  condition  [11], 

(n)  ,  =  |tt(r(().<)  + WO.  1)1/2. 
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The  entropy  condition  may  be  written  as 


(12)  M(t)  =  u,(r(t),  t)  -  t)>0 

for  t  > 

The  outer-region  subdomainB  are  dictated  by  the  regions  where  [/  is  a  good  ap¬ 
proximation  to  u.  These  are  defined  by  bounding  the  difference  U  —  u.  The  boimds 
are  reflected  in  the  following  theorem. 

Theorem  l.  (Howes  (6)).  Let  u{x,t,t)  be  the  solution  to  Pfu]  =  0  on  D  and 
U{x,t)  be  the  solution  to  Po[f^]  =  0  m  the  limit  as  e  tends  to  zero,  each  satisfying  the 
the  boundary  data  (3-5).  Assume  that  the  boundary  data  (3-5)  satisfy  the  compatibility 
eonditiorts  (6), (8-0),  and  that  a,  P  and  with  their  first  and  second  derivatives  are 
all  bounded.  Then  for  e  smalt  enough 

(13)  |tt  -  Ifl  =  0(#iexp[-/*(x,t)/c‘/*|)  +  0{e) 
when  the  denvattves  of  U  are  continuous  across  F,  and 

(14)  |u  -  =  0(/xexp(-/*/c*^*))  +  0{e^^^S  exp[-//c*^*))  -I-  0{e) 

in  the  more  general  ease  when  the  derivatives  ofU  are  not  continuous  across  F.  Here 
f(x,t)  is  a  distance  function  between  (x,t)  and  (F,t),  and  6  is  an  upper  bound  on  the 
difference  of  the  normal  derivative  of  U  across  F. 

The  subdomains  are  dictated  by  the  error  bounds  of  this  theorem.  These  bounds 
are  small  except  in  an  asymptotically  small  neighborhood  of  the  shock.  The  outer 
region  subdomain  is  the  portion  of  D  where  where  using  U  to  approximate  u  intro¬ 
duces  a  small  error.  The  subdomain  where  U  may  be  a  poor  approximation  to  u 
includes  the  internal  or  shock  layer.  The  internal-layer  subdomain  is  the  following 
neighborhood  of  F: 

(15)  Du  =  {(x,t)l(x,t)  eD,\x-  f-MOI  <  A(0>. 

Here  ^(t)  is  the  width  of  the  internal-layer  subdomain  at  time  t.  Theorem  1  dictates 
that  |tt— 17|  =  0{u)  in  D/i  when  the  internal-layer  subdomain  has  size  0[ri[t)i^l*  In^''*  i/). 
Thus,  to  obtain  an  a  priori  bound  of  0(c)  on  the  error,  the  internal  layer  will  be  no 
larger  than  A{t)  <  Krt{t)t^l^h}l^ t,  where  iC  is  a  constant  independent  of  e.  The 
outer-region  subdomain  is  the  complement  of  Du  with  respect  to  D, 

(16)  Doji  =  {(x,*)|(x,t)  €D,ix-  F-*(t)|  >  A(0>- 

Since  the  method  is  designed  for  small  e,  the  internal-layer  subdomun  will  be  an 
asymptotically  small  region  surrounding  the  shock. 

Since  the  solution  to  the  reduced  equation  in  the  shock-layer  subdomain  may 
result  in  large  errors,  the  reduced  equation  will  be  solved  only  in  Dor>  The  method 
will  use  the  full  equation  in  Du  subject  to  data  provided  by  the  solution  to  the  reduced 


equation  in  the  outer  region.  An  analysis  of  the  error  induced  with  this  procedure  is 
presented  in  Corollary  2  below. 

The  local  error  botinds  of  Theorem  1  are  now  used  to  establish  a  global  a  priori 
error  bound  when  using  this  procedure.  The  bound,  as  presented  is  sharp  in  Dqr', 
however,  the  bound  reflects  the  crude  error  bound  of  Theorem  1  in  the  region  of  the 
shock. 

Corollary  2.  Let  u  be  the  solution  to  (l)  satisfying  (3-5).  Suppose  v  is  obtained 
by  first  solving  (7)  m  Dqr  subject  to  (3-S),  then  solving  (l)  on  Dn,  with  boundary  data 
V  on  dDjL.  Assume  that  the  boundary  data  (3-5)  satisfy  the  compatibility  eonditioru 
(6)  ,(8-9),  and  that  a,  0  and  7,  with  their  first  and  second  derivatives  are  all  bounded. 
If  E  =  ||tt  —  v||t,  then  for  e  small  enough 

(17)  E  =  0(c) 


in  Dor,  and 

(18)  E  =  0{e^f*\n^f*e) 


in  Djl. 

In  the  proof,  a  simple  bound  on  the  size  of  the  solution  in  the  internal  layer  is 
established  via  a  maximum  principle  argument.  From  there,  the  proof  follows  directly 
from  Theorem  1. 

Proof.  The  Li  norm  is  deflned  as 

lk(*i*)lli  =  \gix,t)\dxdt. 

Inequality  (17)  follows  directly  from  applying  the  Li  norm  to  the  boimd  (13)  of 
Theorem  1. 

Let 


w(t)  =  K\  +  K^t. 

The  constant  Kt  may  be  chosen  independent  of  c  such  that 

P[w]  >  P[w]  =  0 

for  (x,()  €  Dii,  and  the  constant  Ki  may  be  chosen  independent  of  e  such  that 

w  >  w 

for  (x,  t)  €  BDii..  Thus,  the  conditions  of  the  statement  of  the  maximum  principle 
due  to  Nagumo  and  Westfal  [18]  are  satisfied,  and  u;  >  v  for  all  of  Djt.  A  symmetric 
argument  can  be  used  to  establish  a  lower  boimding  function  which  has  the  same 
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fonn  as  w.  In  addition,  a  similar  argument  can  be  used  to  establish  upper  and  lower 
bounds  for  u  in  D/f  Thus, 


(tt  -  v]  =  0(1) 

for  (x,t)  6  Dji.  Since  the  area  of  D/i  is  0(c*/*  In'^*  f),  applying  the  Li  norm,  the 
bound  (18)  follows.  □ 

Thus  far,  the  problem  has  been  preconditioned  by  decomposing  the  domain  into 
regions  where  the  solution  behaves  differently,  and  forming  subproblems  in  those 
re^ons.  Another  preconditioning  for  the  problem  in  the  internal-layer  subdomain  is 
also  appropriate. 

The  preconditioning  in  the  internal  layer  is  a  scaling  and  translation  of  the  spatial 
coordinate.  The  translation  allows  the  coordinate  system  to  move  with  the  shock  by 
using  X  —  r(t)  in  place  of  x.  The  scaling  is  to  stretch  the  spatial  coordinate  by  l/e. 
This  scaling  is  identified  via  multiple-scale  asymptotic  analysis  [9,12,16],  and  allows 
the  shock  to  be  resolved  in  the  local  coordinate  systexn.  Combining  these  two,  the 
new  spatial  coordinate  in  Du  is 

(19)  x=(x  — r)/c. 

The  computational  analog  of  this  transformation  is  described  in  [15]  or  [10]. 

Asymptotics  identified  two  subdomains  and  provided  preconditioners  for  the  prob¬ 
lems  within  both  subdomains.  The  domain  decomposition  algorithm  described  thus 
far  requires  a  priori  knowledge  of  the  location  of  the  shock.  To  remove  this  restriction, 
the  domain  decomposition  is  combined  with  a  functional  iteration.  This  allows  for 
the  iterative  determination  of  the  shock  location  in  the  computational  method. 

4.  Iteration.  The  domun  decomposition  algorithm  described  in  the  previous 
section  will  now  be  treated  as  a  single  step  in  an  iterative  process.  Each  step  of 
this  iteration  requires  the  solution  of  a  linearized  form  of  the  reduced  equation  (7)  in 
the  outer-region  subdomain  followed  by  the  solution  of  the  full  equation  (1)  in  the 
internal-layer  subdomain.  Denote  the  iterate  by  The  equations  governing  the 

iterate  are 

(20)  =  0 

in  Dor,  and 

(21)  0^*^  +  =  0 

in  Dji,  Boundary  data  for  the  internal-layer  subdomain  is  provided  by  the  solution 
of  (20)  in  the  outer-repon  subdomain. 

In  this  section,  the  convergence  of  the  iteration  (20)  in  the  outer-region  subdomain 
to  a  solution  of  (7)  will  be  established.  In  addition,  a  global  o  priori  error  bound  for 
the  method  will  be  presented.  Throughout  this  section  the  conditions  on  the  boundary 
data  presented  in  Section  2  are  assumed  to  be  satisfied. 

S 


Theorem  3.  let  be  the  set  of  iterates  of  (20)  m  the  subdomain 

Dor  satisfying  the  boundary  data  (3-5)  with  initial  guess  Assume  satisfies 
(3-5)  and  t«  Lipsehitz  continuous  on  D.  Let 

S  =  sup  —  ^*'*1. 

D 

Then 

(22)  il7*+^  -  t;*l  <  SCe'^*{e^  -  1) 

for  (z,  t)  €  Dor,  where  C,  X  and  R  are  known  positive  constants. 

The  proof  utilizes  some  results  on  continuity  of  the  iterates  which  will  be  estab¬ 
lished  first.  The  boundedness  of  is  the  subject  of  the  following  lemma. 

Lemma  4.  £et  be  the  solution  to  equation  (20)  on  the  subdomain  Dqr, 
where  is  Lipsehitz  continuous.  Then, 

<  A’e", 

for  (z,  t)  €  Dor,  where  K  and  R  are  constants  independent  of  x,  t  and  e. 

Proof.  Consider  the  transformation  (z,t)  — » ({,r)  defined  by 

(23)  t  =  r, 
and 

(24)  = 
with  initial  conditions 

(25)  x‘(0)  =  e.  5>e>0; 

(26)  **(yo  MO)  =  0,  €  <  0; 

(27)  z*(y»-M0)  =  0,  i>b. 

Here,  (^,r)  =  (yo(»’)f  0  **  image  of  the  curve  (z,t)  =  (0,t),  and  ({,r)  =  (y»(r),T) 
is  the  image  of  the  curve  (z,t)  =  (6,t).  Under  this  transformation,  equation  (20) 
becomes 

(28) 

Since  is  Lipsehitz  continuous,  the  transformation  (23)-(27)  is  uniquely  defined, 
and  equation  (28)  may  be  used  in  place  of  equation  (20). 
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The  coefficient  r  is  bounded,  hence  there  is  a  constant  fC  independent  of  z,  t,  and 
€  such  that  O  =  A’e*'  is  an  upper  bound  for  In  addition,  ^  =  —O  is  a  lower 

bound.  Since  t  =  r  the  desired  result  is  established.  □ 

Based  on  the  assumption  of  the  boundedness  of  the  boundedness  oi  is 
established  in  Lemma  5. 

Lemma  5.  Suppott  that  the  eonditioru  of  Lemma  4  obtain.  Then  0^*^  and 
are  bounded  independent  of  x,  t  and  e  in  Dqr- 

In  the  proof  of  this  lemma,  an  equation  governing  will  be  derived.  Then,  a 
form  of  the  maximum  principle  will  be  shown  to  apply  after  a  change  of  the  dependent 
variable. 

Proof.  Let  w  —  0^*'.  The  equation 

(29)  w,  =  (r  -  Ui)w  + 

is  derived  by  taking  the  partial  with  respect  to  z  of  equation  (20),  then  applying  the 
transformation  (23)-(27).  Boundary  data  for  w  may  be  obtained  by  differentiating 
(3-5). 

Let  w  —  e^'v.  The  equation  governing  v  is 

t;,  =  (r-lf.‘-A)v+^t/‘. 

ax 

Choose  A  =  -  max(0,  inf /)(r(z)  -  Ojl)),  so  that  the  coefficient  of  v  in  this  equation 
will  be  nonnegative. 

Define  an  upper  bounding  function  as 

where  the  constants  Ki  and  Kt  will  be  chosen.  Then  z  =  0  -  v  satisfies 

(30)  Xr  =  A(x,t)x  +  B(z,f), 

where  A  =  r  —  0^  —  X  and  B  —  [Ifj  —  {r  —  O^  —  A)]C»  +  The  constants 

Ki  and  may  be  chosen  so  that  A,  B  and  z(z,  0)  are  nonnegative.  For  example, 
chooee 

Kj  =  max(2,sup(r  -  -  A))  and  Ki  =  max(0,sup  ^&*,sup2w). 

Don  Don  »*  n 

Under  these  conditions,  z  is  poeitive,  and  O  is  an  upper  bound  for  v. 

A  lower  bounding  function  may  be  obtuned  in  a  nmilar  way.  Set  u.=  Ki  (e^*^  -  | 
If  Ki  and  A  are  chosen  as  before  and 

Kt  =  min(0,  inf  ^^*,inf  2w), 

Don  dx  n 
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then  —  V  is  nonpoeitive.  Thus,  is  a  lower  bounding  function  for  v. 

Since  r  <  T,  the  boundedneae  of  independent  of  x,  t  and  e  follows.  The 
boundedness  of  follows  from  the  boundedness  of  the  terms  in  equation  (20). 
Therefore,  the  desired  result  is  established.  □ 

Given  Jiat  is  Lipschits  continuous,  Lemma  5  states  that  all  of  the  remaining 
interates  will  be  Lipschits  continuous.  The  significance  of  this  is  that  the  characteristic 
transformation  (23)-(27)  is  uniquely  defined.  This  is  used  in  the  proof  of  the  theorem. 

Proof.  (Theorem  3).  The  equation  governing  *  =  t/*  -  y*'*’*  in  the  characteristic 
coordinate  system  (23-27)  is 


+  rz. 

An  exponential  change  of  variable  will  transform  the  problem  such  that  a  form 
of  the  manmum  principle  applies.  Let 

z  = 


where  A  =  -  max(0,inf  r(x)).  Thus,  the  equation  for  u>  is 

w,  =  «*‘(t/*-‘  -  +  ho. 


where  ^  =  r  +  A,  and  w  =  0  on  11. 

An  upper  bound  for  w  may  now  be  defined.  Let 


u;(r)  =  6 


n 

/f(e*t-l)’ 


where  R  =  supr  -  inf  r  =  supf.  The  constant  q  will  be  determined  shortly.  Taking 
the  partial  derivative  with  respect  to  r. 


w,  =  tfq  +  R(jj. 


An  equation  for  /  =  w  —  tv  is 

(31)  /,  =  \6ri  -  e*'(&*-‘  -  -»-(«-  f)^]  f /. 

Since  w  >  0  and  tv  =  0  on  11,  the  fimction  /  is  nonnegative  on  IT.  Choose  q  = 
max(0,  info  Thus,  the  source  term  in  equation  (31)  is  nonnegative,  and 

/  >  0  in  DoM’  Thus  w(f)  >  tv((,r)  for  (x({,r),l((,r))  in  Dor-  But  t  =  t,  hence 


d#fe-**(e**  -  1)  >  s  =  -  0^, 


for  K  —  ti/R. 

Using  symmetric  arguments,  —u  can  be  shown  to  be  a  lower  bound  on  z,  leading 
to  the  desired  result.  □ 
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This  theorem  provides  an  upper  bound  on  the  latest  time  for  which  the  iteration 
converges.  Apply  the  infinity  norm  to  (22)  to  obtain 

Then  the  following  corollary  provides  the  conditions  for  convergence. 

Corollary  6.  Suppose  that  the  eonditiona  of  Theorem  3  obtain.  Reatriet  the 
domain  to  Don-  be  the  largest  positive  number  such  that 

C  =  sup  Ce-*‘(e*‘  -  1)  <  1. 

Then  the  sequence  of  iterates  defined  by  (20)  converges  to  a  solution  of  (7)  satisfying 
the  data  (3-5)  for  the  upper  bound  T  on  time  of  equation  (2)  satisfying  Tmnx  >  T  >  0. 

Proof.  With  the  restriction  of  t  <  Taaxt  the  iteration  is  a  contraction  mapping, 
and  the  result  follows.  □ 

A  statement  of  a  global  a  priori  error  bound  for  the  computational  method  is 
presented  in  Corollary  7  below.  As  with  Corollary  2,  the  bound  is  sharp  in  Don', 
however,  the  bound  is  crude  in  the  region  of  the  shock. 

Corollary  7.  Let  u  be  the  solution  to  (l)  satisfying  (3-5).  Suppose  each  iterate 
17*  obtained  by  first  solving  (20)  m  Dor  subject  to  (3-5),  then  solving  (21)  on  Du 
vnth  boundary  data  on  dDji.  Suppose  T^  >  T  >  0,  and  let  v  =  limt_oo  = 
O’**.  If  E  =  Iju  -  vlli,  then  for  e  small  enough 

(32)  E  =  0{e) 
in  Dor,  and 

(33)  JS;  =  0(€‘/*ln*/*«) 
m  Dji.  Here  A  =  0(c‘/Mn‘/’  «)• 

To  prove  Corollary  7,  it  will  first  be  established  that  O’**  is  the  desired  solution 
of  equation  (7).  The  bound  in  (32)  will  follow.  Then,  Lemma  8  will  be  applied  to 
establish  the  bound  (33). 

LEMMA  8.  Let  v  be  the  solution  to  equation  (1)  on  Du,.  Suppose  that  the  data 
specified  on  dD/i  are  bounded  with  bounded  derivatives.  Assume  that  dDn  is  at  least 
C>.  Then 

IH  <  K, 

for  (z,<)  e  Dil, 

In  the  proof  of  this  lemma,  upper  and  lower  bounding  functions  for  v  are  estab¬ 
lished  by  applying  the  version  of  the  maximum  principle  due  to  Nagumo  and  Westfal 
118). 
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Proof.  Make  the  change  of  dependent  variable  w  =  where  X  satisfies 

r  =  r  —  A<fo<0.  Then 


:=  Wt  +  -  rw  —  £w„  =  0. 

where  w  =  e~**v  on  dDn.  Define 


Then 


G)  =  k  =  xnax(0,  sup 

SDii, 


Px[0\  =  -foAT  >  PxlH. 

In  addition,  u»  >  e~*‘*v  =  to  for  (i,  t)  €  dDjj,.  These  conditions  allow  the  application 
of  the  Nagumo-Westfal  Lemma  to  conclude  Q  >  w  for  (®,  t)  €  Du.  By  symmetric 
arguments,  jit  =  —Q  may  be  shown  to  be  an  lower  bounding  function  on  to.  Setting 
K  =  e^^kf  the  result  is  established.  □ 

Proof  (Corollary  7),  In  the  outer  region,  v  =  where  bounds  for  u  —  w  = 
u  —  are  given  in  (13-14).  Thus,  there  is  a  constant  ffoA  such  that  |u- 1/^*1  <  eKoR 
for  (x,  t)  €  Dor.  Applying  the  Li  norm  to  u  —  O'”,  relation  (32)  is  established. 

From  Lemma  8,  there  is  a  constant  Kn,  such  that  |u|  +  |o|  <  Kn  in  Du,  where 
K;i  is  independent  of  e.  Thus, 


Since  the  area  covered  by  Dil  is  of  size  0{e^^*  In}^*  e) ,  relation  (33)  holds.  □ 

5.  Concluding  Remarks.  A  computational  method  will  be  created  from  the 
theory  presented  in  the  past  few  sections.  The  method  is  independent  of  the  par¬ 
ticular  numerical  schemes  used;  however,  candidates  for  the  numerical  schemes  will 
be  discussed.  The  computational  method  is  be  constructed  by  solving  equation  (21) 
then  solving  equation  (20)  in  succession.  With  minor  modifications  to  account  for  the 
forcing  term,  the  method  of  characteristics  may  be  used  for  (21),  or  the  equation  my 
be  solved  using  one  of  methods  discussed  in  [1,17,3,14].  Equation  (20)  is  solved  in  the 
locsd  coordinate  system  (10).  In  this  coordinate  system,  the  coefficient  of  the  diffusion 
is  large  enough  to  allow  the  application  of  standard  finite  difference  methods. 

In  the  most  general  case,  the  internal  layer  subdomain  will  change  from  iteration 
to  iteration.  The  theory  presented  in  tlus  paper  restricts  consideration  to  a  stationary 
internal-layer  subdomain;  however,  experimental  results  demonstrate  that  this  is  not 
a  constr^nt  in  the  numerical  method  [lOj.  In  addition,  extensions  to  the  theoretical 
basis  to  include  a  moving  boundary  are  the  subject  of  pronusing  current  research. 
When  dDji  is  allowed  to  move  between  iterations,  the  method  must  be  able  to  deter¬ 
mine  the  location  of  the  boundary  as  the  computations  proceed.  This  may  be  done 
by  monitoring  the  Jacobian  of  the  characteristic  transformation  (23-27),  as  discussed 
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in  [10].  In  addition,  Lemma  5  suggests  that  monitoring  may  be  used  to  place  the 
boundary.  The  value  of  will  be  large  and  negative  in  a  neighborhood  of  the  shock, 
and  it  will  be  bounded  independent  of  6  in  the  outer-region. 

This  paper  demonstrated  the  use  of  asymptotics  to  dictate  a  numerical  method 
with  high  accuracy  and  efficiency.  Asymptotic  analysis  provided  a  theoretical  basis 
for  a  dom^  decomposition,  and  guided  in  the  derivation  of  rigorous  local  and  global 
a  priori  error  bounds. 

The  asymptotic  analysis  can  be  used  to  analyze  existing  methods  as  well  as  de¬ 
velop  new  algorithms.  For  example,  the  analysis  may  be  used  to  study  the  effects  of 
using  artificial  viscosity.  Theorem  1  provides  an  a  priori  upper  bound  on  the  error 
induced  by  using  i  in  place  of  e  where  2  >  c.  Such  an  equation  would  be  solved  when 
uring  constant  coefficient  artificial  diffusion.  This  bound  may  be  obtained  by  using 
the  theorem  to  establish  bounds  separately  on  |u(z,  t,  e)  —  U\  and  |u(x,  t,  e)  —  U\,  then 
summing  them. 

The  availability  of  estimates  and  bounds  on  the  error  is  important  in  the  design  of 
numerical  methods.  Rigorous  a  priori  error  bounds  were  established  for  the  method 
presented  here.  In  addition,  the  particular  numerical  schemes  used  for  the  subprob¬ 
lems  allowed  a  posfertbrt  error  estimation.  The  a  prton  error  bounds  were  shown 
to  be  much  larger  than  the  errors  observed  in  the  computations;  thus,  sharper  error 
botmds  are  expected. 
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